// Copyright 2020 Poofee (https://github.com/Poofee)
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// ------------------------------------------------------------------------
/*****************************************************************************
 *                                                                           *
 *  File:    relay.cpp                                                       *
 *                                                                           *
 *****************************************************************************
 *                                                                           *
 *  Authors: Poofee                                                          *
 *  Email:   poofee@qq.com                                                   *
 *  Address: https://github.com/Poofee                                       *
 *  Date:    2020年07月15日                                                   *
 *                                                                           *
 *****************************************************************************/
#include "relay.h"


#include <time.h>

Relay::Relay() :
    current_step(0),
    ParallelThreads(1),
    unitRatio(1e-3),    
	pmeshnode(nullptr),
	pmeshele(nullptr),
	pmeshnodelast(nullptr),
	pmeshelelast(nullptr),
    model(nullptr),
	RemeshOrMorph(true) {
    setPrefixName((char*)"");
}

Relay::~Relay() {
	if(model) gmsh::finalize();

    if (pmeshele)      {delete pmeshele;      pmeshele=nullptr;}
    if (pmeshnode)     {delete pmeshnode;     pmeshnode=nullptr;}
    if (pmeshelelast)  {delete pmeshelelast;  pmeshelelast=nullptr;}
    if (pmeshnodelast) {delete pmeshnodelast; pmeshnodelast=nullptr;}

	for (auto m : materialList) {
        if (m) {
            delete m;
            m = nullptr;
        }
	}
	materialList.clear();
	materialMap.clear();
}


/*!
 * @brief      初始化函数。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     void  
**/
void Relay::init() {

}


/*!
 * @brief      设置文件名。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  char * fn  
 * @return     void  
**/
void Relay::setFileName(char* fn) {
	sprintf(fileName, "%s", fn);
}


/*!
 * @brief      设置单位到m。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  double ratio  
 * @return     void  
**/
void Relay::unitConvert(double ratio) {
	for (int i = 0; i < num_pts; ++i) {
		pmeshnode[i].x *= ratio;
		pmeshnode[i].y *= ratio;
		pmeshnode[i].z *= ratio;
	}
}


/*!
 * @brief      打开geo文件。注意，由于瞬态问题需要使用到上一次的
               数据，所以，最好在geo文件中，将发生形变的区域移到最后。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     void  
**/
void Relay::openGeo() {
    printCurrentTime();
	clock_t time;
	time = SuperLU_timer_();

	int myargn = 4;
    char geoName[256];
    sprintf(geoName, "%s.geo", fileName);
	char *myargv[] = { (char*)"gmsh",(char*)"-format",(char*)"msh2",(char*)"-v",(char*)"0" };
	gmsh::initialize(myargn, myargv);
	gmsh::option::setNumber("General.Terminal", 1);
	gmsh::open(geoName);
	//    gmsh::model::mesh::generate(2);
	//    GModel::current()->mesh(2);
	model = GModel::current();

    printf("Opening model %s\n\n", geoName);
    printCurrentTime();
	printf("Meshing...\n");
	/** 初次分网，重分网必须先有一个网格 **/
	char nameMsh[256];
	GenerateMesh(model, 3);
	sprintf(nameMsh, "%s_%02d.msh", fileName, 0);
	gmsh::write(nameMsh);

	printf("Mesh time elasped %lf s\n", SuperLU_timer_() - time);
}
int next_int(char **start) {
	int i;
	char *end;

	i = strtol(*start, &end, 10);
	*start = end;
	return(i);
}

/*!
 * @brief      读取当前时间步的mesh，如果当前含有分网数据，先删除。
 gmsh网格的版本为22。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     void  
**/
void Relay::loadMesh() {
	/** 先删除上一次的网格 **/
	if (pmeshelelast) delete pmeshelelast;
	if (pmeshnodelast) delete pmeshnodelast;
	/** 保存上一次的网格 **/
	if (pmeshele) pmeshelelast = pmeshele;
	if (pmeshnode) pmeshnodelast = pmeshnode;

	char nameMsh[256];
	sprintf(nameMsh, "%s_%02d.msh", fileName, current_step);
	char *ch = (char *)calloc(256, sizeof(char));
	//------------open file----------------------------------
	FILE * fp = nullptr;
	fp = fopen(nameMsh, "r");
	if (fp == nullptr) {
		printf("Error: openning file %s!", nameMsh);
        exit(-1);
	}
	while (!feof(fp)) {
		fgets(ch, 256, fp);

		if (strstr(ch, "$MeshFormat")) {
			double version;
			int file_type;
			int data_size;
			if (fscanf(fp, "%lf %d %d\n", &version, &file_type, &data_size) != 3) {
				printf("error reading format");
				return;
			} else {
				if (version > 2.2) {
					printf("Can only open gmsh version 2.2 format");
					return;
				}
			}
			fgets(ch, 256, fp);
			if (!strstr(ch, "$EndMeshFormat")) {
				printf("$MeshFormat section should end to string $EndMeshFormat:\n%s\n", ch);
			}
		} else if (strstr(ch, "$Nodes")) {
			if (fscanf(fp, "%d\n", &num_pts) != 1) {
				return;
			} else {
				/** 读取节点坐标 **/
				pmeshnode = new CNode[num_pts];
				int index;
				for (int i = 0; i < num_pts; ++i) {
					fscanf(fp, "%d %lf %lf %lf\n", &index, &pmeshnode[i].x, &pmeshnode[i].y, &pmeshnode[i].z);
					//qDebug()<<index<<pmeshnode[i].x<<pmeshnode[i].y<<pmeshnode[i].z;
				}
			}
			fgets(ch, 256, fp);
			if (!strstr(ch, "$EndNodes")) {
				printf("$Node section should end to string $EndNodes:\n%s\n", ch);
			}
		} else if (strstr(ch, "$Elements")) {
			int ele_number;
			//    int elm_type;
			int number_of_tags;
			char * chtmp;
			if (fscanf(fp, "%d\n", &num_ele) != 1) {
				return;
			} else {
				pmeshele = new CElement[num_ele];
				for (int i = 0; i < num_ele; ++i) {
					chtmp = fgets(ch, 256, fp);
					ele_number = next_int(&chtmp);
					pmeshele[i].ele_type = next_int(&chtmp);
					number_of_tags = next_int(&chtmp);
					pmeshele[i].physic_tag = next_int(&chtmp);
					pmeshele[i].geometry_tag = next_int(&chtmp);
					pmeshele[i].miu = miu0;
					pmeshele[i].miut = miu0;

					int element_pmeshnode = 0;
					switch (pmeshele[i].ele_type) {
					case 15:/** point **/
						element_pmeshnode = 1;
						break;
					case 1:/** line **/
						element_pmeshnode = 2;
						break;
					case 2:/** triangle **/
						element_pmeshnode = 3;
						break;
					case 3:/** quadrangle **/
						element_pmeshnode = 4;
						break;
					case 4:/** tetrahedron **/
						element_pmeshnode = 4;
						break;
					default:
						element_pmeshnode = 0;
						break;
					}

					for (int j = 0; j < element_pmeshnode; ++j)
						pmeshele[i].n[j] = next_int(&chtmp) - 1;
				}
			}
			fgets(ch, 256, fp);
			if (!strstr(ch, "$EndElements")) {
				printf("$Element section should end to string $EndElements:\n%s\n", ch);
			}
		}
	}
	delete ch;
	fclose(fp);
}


/*!
 * @brief      默认二维查找外部边界，确定dof。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int index  
 * @return     void  
**/
void Relay::findBoundaryPoints2D(int index) {
	/** 根据边界生成所有的点 **/
	int num_bndpoints = boundaryEdges.size() * 2;
	boundaryPoints.resize(num_bndpoints);
	for (int i = 0; i < boundaryEdges.size(); ++i) {
		boundaryPoints[i * 2 + 0] = boundaryEdges[i].start;
		boundaryPoints[i * 2 + 1] = boundaryEdges[i].end;
	}
	/** 去重 **/
	std::vector<int>::iterator it_1 = boundaryPoints.begin();
	std::vector<int>::iterator it_2 = boundaryPoints.end();
	std::vector<int>::iterator new_end;

	sort(it_1, it_2);
	new_end = unique(it_1, it_2);
	boundaryPoints.erase(new_end, it_2);
	printf("Total %llu boundary points.\n", boundaryPoints.size());
}


/*!
 * @brief      查找二维模型外边界上的边
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int index  查找区域index，如果未指定，则是全局。
 * @return     void  
**/
void Relay::findBoundaryEdges2D(int index) {
	/** 计算出三角单元的数目 **/
	num_triangle = 0;
	for (int i = 0; i < num_ele; ++i) {
		if (pmeshele[i].ele_type == 2) {
			num_triangle++;
		}
	}
	/** 生成所有的棱 **/
	boundaryEdges.resize(num_triangle * 3);
	allPoints.resize(num_triangle * 3);

	for (int i = 0; i < num_triangle; ++i) {
		int a, b, c, t;
		int i1 = num_ele - num_triangle + i;
		a = pmeshele[i1].n[0];
		b = pmeshele[i1].n[1];
		c = pmeshele[i1].n[2];
		/** 按照升序排序 **/
		if (a > b) { t = a; a = b; b = t; }
		if (a > c) { t = a; a = c; c = t; }
		if (b > c) { t = b; b = c; c = t; }

		boundaryEdges[i * 3 + 0].start = a;
		boundaryEdges[i * 3 + 0].end   = b;
		boundaryEdges[i * 3 + 0].bdr   = 1;
		boundaryEdges[i * 3 + 1].start = a;
		boundaryEdges[i * 3 + 1].end   = c;
		boundaryEdges[i * 3 + 1].bdr   = 1;
		boundaryEdges[i * 3 + 2].start = b;
		boundaryEdges[i * 3 + 2].end   = c;
		boundaryEdges[i * 3 + 2].bdr   = 1;

		allPoints[i * 3 + 0] = a;
		allPoints[i * 3 + 1] = b;
		allPoints[i * 3 + 2] = c;
	}

	/** 去重找到边界 **/
	std::vector<FEMedge>::iterator it_1 = boundaryEdges.begin();
	std::vector<FEMedge>::iterator it_2 = boundaryEdges.end();
	std::vector<FEMedge>::iterator new_end;

	std::sort(it_1, it_2);

	new_end = Myunique(it_1, it_2);
	allEdges = boundaryEdges;
	boundaryEdges.erase(new_end, it_2);
	//    for(int k = 0;k < boundaryEdges.size();++k){
	//        printf("%d: start: %d,end:%d\n",k,boundaryEdges[k].start,boundaryEdges[k].end);
	//    }
	/** 去重找到所有的节点 **/
	std::vector<int>::iterator it_3 = allPoints.begin();
	std::vector<int>::iterator it_4 = allPoints.end();
	std::vector<int>::iterator new_end1;

	std::sort(it_3, it_4);

	new_end1 = std::unique(it_3, it_4);
	allPoints.erase(new_end1, it_4);

	printf("Total %d triangles.\n", num_triangle);
	printf("Total %llu points.\n", allPoints.size());
	printf("Total %llu boundary edges.\n", boundaryEdges.size());
}

/*!
 * @brief      三维版的查找边界棱。计算出三维分网棱的连接关系。
               三维的话，要计算边界棱的话，还得先计算出边界面，感觉还是
               应该计算出连接关系。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int index  
 * @return     void  
**/
void Relay::findBoundaryEdges3D(int index) {
    double time;
	time = SuperLU_timer_();
	
    /** 去重找到边界 **/
	std::vector<FEMedge>::iterator it_1 = allEdges.begin();
	std::vector<FEMedge>::iterator it_2 = allEdges.end();
	std::vector<FEMedge>::iterator new_end;

	std::sort(it_1, it_2);

	new_end = unique(it_1, it_2);
	//allEdges = boundaryEdges;
	allEdges.erase(new_end, it_2);
	/** 重分网后hashmap需不需要初始化？看起来貌似不需要，因为直接就替换了 **/
	size_t len = allEdges.size();
	for (size_t i = 0; i < len; i++) {
		FEMedge e = allEdges[i];
		alledgeshash[std::pair<int,int>(e.start, e.end)] = i;
	}

	/** 建立连通矩阵 **/
	int n[4];
    for (int i = 0; i < num_tet; ++i) {
        CElement* ele = pmeshele + num_ele - num_tet + i;

		for (int in = 0; in < 4; ++in) {
			n[in] = ele->n[in];
		}
		int start, end;
		for (int in = 4; in < 10; ++in) {
			start = n[tet_n[in-4][0]];
			end   = n[tet_n[in-4][1]];
			std::unordered_map<std::pair<int, int>, int, pair_hash>::iterator it;
			if (start < end) {
				it = alledgeshash.find(std::pair<int, int>(start, end));
			} else {
				it = alledgeshash.find(std::pair<int, int>(end, start));
			}
			
			if (it != alledgeshash.end()) {
				ele->n[in] = it->second;
			} else {
				printf("Error in void Relay::findBoundaryEdges3D(int index)\n");
                exit(-1);
			}
		}
	}
	/** 查找边界棱 **/
	int bndFaceSize = boundaryFaces.size();
	boundaryEdges.resize(bndFaceSize * 3);
	boundaryPoints.resize(bndFaceSize * 3);
	int a, b, c, t;
	for (int iF = 0; iF < bndFaceSize; ++iF) {
		auto f = boundaryFaces[iF];
		a = f.n1;
		b = f.n2;
		c = f.n3;
		/** 按照升序排序 **/
		if (a > b) { t = a; a = b; b = t; }
		if (a > c) { t = a; a = c; c = t; }
		if (b > c) { t = b; b = c; c = t; }

		boundaryEdges[iF * 3 + 0].start = a;
		boundaryEdges[iF * 3 + 0].end   = b;
		boundaryEdges[iF * 3 + 1].start = a;
		boundaryEdges[iF * 3 + 1].end   = c;
		boundaryEdges[iF * 3 + 2].start = b;
		boundaryEdges[iF * 3 + 2].end   = c;

		boundaryPoints[iF * 3 + 0] = a;
		boundaryPoints[iF * 3 + 1] = b;
		boundaryPoints[iF * 3 + 2] = c;
	}

	/** 去重找到边界 **/
	std::vector<FEMedge>::iterator it_1b = boundaryEdges.begin();
	std::vector<FEMedge>::iterator it_2b = boundaryEdges.end();
	std::vector<FEMedge>::iterator new_endb;

	std::sort(it_1b, it_2b);

	new_endb = unique(it_1b, it_2b);
	boundaryEdges.erase(new_endb, it_2b);

	/** 计算边界点 **/
	std::vector<int>::iterator it_1p = boundaryPoints.begin();
	std::vector<int>::iterator it_2p = boundaryPoints.end();
	std::vector<int>::iterator new_endp;

	sort(it_1p, it_2p);
	new_endp = unique(it_1p, it_2p);
	boundaryPoints.erase(new_endp, it_2p);

    num_edges = allEdges.size();

    printf("Total %d edges.\n", num_edges);
	printf("Total %llu boundary edges.\n", boundaryEdges.size());
	printf("Total %llu boundary points.\n", boundaryPoints.size());
	printf("findBoundaryEdges3D time elasped %lf s\n", SuperLU_timer_()-time);
}

/*!
 * @brief      查找三维体的边界上的面单元。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int index  为-1时，计算整体的边界。
 * @return     void  
**/
void Relay::findBoundaryFaces3D(int index) {
	/** 计算出四面体单元的数目 **/
	num_tet = 0;
	for (int i = 0; i < num_ele; ++i) {
		if (pmeshele[i].ele_type == 4) {
			num_tet++;
		}
	}
	/** 生成所有的面 **/
	boundaryFaces.resize(num_tet * 4);
	allEdges.resize(num_tet * 6);
	allPoints.resize(num_tet * 4);

	for (unsigned int i = 0; i < num_tet; ++i) {
		int a, b, c, d, t;
		int i1 = num_ele - num_tet + i;
		CElement* ele = pmeshele + i1;
		a = ele->n[0];
		b = ele->n[1];
		c = ele->n[2];
		d = ele->n[3];
		/** 按照升序排序 **/
		if (a > b) { t = a; a = b; b = t; }
		if (a > c) { t = a; a = c; c = t; }
		if (a > d) { t = a; a = d; d = t; }
		if (b > c) { t = b; b = c; c = t; }
		if (b > d) { t = b; b = d; d = t; }
		if (c > d) { t = c; c = d; d = t; }

		allEdges[i * 6 + 0].start = a;
		allEdges[i * 6 + 0].end   = b;
		allEdges[i * 6 + 0].bdr   = 0;
		allEdges[i * 6 + 1].start = a;
		allEdges[i * 6 + 1].end   = c;
		allEdges[i * 6 + 1].bdr   = 0;
		allEdges[i * 6 + 2].start = a;
		allEdges[i * 6 + 2].end   = d;
		allEdges[i * 6 + 2].bdr   = 0;
		allEdges[i * 6 + 3].start = b;
		allEdges[i * 6 + 3].end   = c;
		allEdges[i * 6 + 3].bdr   = 0;
		allEdges[i * 6 + 4].start = b;
		allEdges[i * 6 + 4].end   = d;
		allEdges[i * 6 + 4].bdr   = 0;
		allEdges[i * 6 + 5].start = c;
		allEdges[i * 6 + 5].end   = d;
		allEdges[i * 6 + 5].bdr   = 0;

		//a,b,c
		boundaryFaces[i * 4 + 0].n1 = a;
		boundaryFaces[i * 4 + 0].n2 = b;
		boundaryFaces[i * 4 + 0].n3 = c;
		//a,b,d
		boundaryFaces[i * 4 + 1].n1 = a;
		boundaryFaces[i * 4 + 1].n2 = b;
		boundaryFaces[i * 4 + 1].n3 = d;
		//a,c,d
		boundaryFaces[i * 4 + 2].n1 = a;
		boundaryFaces[i * 4 + 2].n2 = c;
		boundaryFaces[i * 4 + 2].n3 = d;
		//b,c,d
		boundaryFaces[i * 4 + 3].n1 = b;
		boundaryFaces[i * 4 + 3].n2 = c;
		boundaryFaces[i * 4 + 3].n3 = d;

		allPoints[i * 4 + 0] = a;
		allPoints[i * 4 + 1] = b;
		allPoints[i * 4 + 2] = c;
		allPoints[i * 4 + 3] = d;
	}

	/** 去重找到边界 **/
	std::vector<FEMface>::iterator it_1 = boundaryFaces.begin();
	std::vector<FEMface>::iterator it_2 = boundaryFaces.end();
	std::vector<FEMface>::iterator new_end;

	std::sort(it_1, it_2);

	new_end = Myunique(it_1, it_2);
	allFaces = boundaryFaces;
	boundaryFaces.erase(new_end, it_2);

	/** 去重找到所有的节点 **/
	std::vector<int>::iterator it_3 = allPoints.begin();
	std::vector<int>::iterator it_4 = allPoints.end();
	std::vector<int>::iterator new_end1;

	std::sort(it_3, it_4);

	new_end1 = std::unique(it_3, it_4);
	allPoints.erase(new_end1, it_4);

	printf("Total %d tetrahedrals.\n", num_tet);
	printf("Total %llu points.\n", allPoints.size());
	printf("Total %llu boundary faces.\n", boundaryFaces.size());
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int xiantie  
 * @return     void  
**/
void Relay::setXiantieTag(int xiantie) {
	tag_xiantie = xiantie;
}


/*!
 * @brief      用来最后输出vtk的时候，不显示空气区域的网格。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int air  
 * @return     void  
**/
void Relay::setAirTag(int air) {
	tag_air = air;
}


/*!
* @brief      
* @author     Poofee
* @date       2020年8月25日
* @param[out] 
* @param[in]  int remesh  
* @return     void  
**/
void Relay::setRemeshTag(int remesh) {
	tag_remesh = remesh;
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  GFace * f  
 * @return     void  
**/
void Relay::deleteFaceMesh(GFace *f) {
	deMeshGFace dem;
	dem(f);
}


/*!
 * @brief      删除区域内的网格。应该是内部的网格。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  GRegion * r  
 * @return     void  
**/
void Relay::deleteRegionMesh(GRegion* r) {
	deMeshGRegion dem;
	dem(r);
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int timestep  
 * @return     void  
**/
void Relay::remesh(int timestep) {
	printf("Must be implemented in subclass.\n");
}


/*!
 * @brief      位移为正值，沿着轴的正方向运动，反之，为负。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  GFace * f  
 * @param[in]  double dx  
 * @param[in]  double dy  
 * @param[in]  double dz  
 * @return     void  
**/
void Relay::moveFace(GFace *f, double dx, double dy, double dz) {
	/** 修改内部分网节点坐标 **/
	for (std::size_t j = 0; j < f->getNumMeshVertices(); ++j) {
		f->getMeshVertex(j)->setXYZ(f->getMeshVertex(j)->x() + dx,
			f->getMeshVertex(j)->y() + dy,
			f->getMeshVertex(j)->z() + dz);
	}
	/** 修改边界棱上分网节点坐标 **/
	std::set<MVertex *, MVertexLessThanNum> all_vertices;
	for (auto e : f->edges()) {
		for (auto line : e->lines) {
			MVertex *v1 = line->getVertex(0);
			MVertex *v2 = line->getVertex(1);

			all_vertices.insert(v1);
			all_vertices.insert(v2);
		}
	}
	/** all_vertices比e->getMeshVertex多了顶点 **/
	for (std::set<MVertex *, MVertexLessThanNum>::iterator ite = all_vertices.begin(); ite != all_vertices.end(); ite++) {
		(*ite)->setXYZ((*ite)->x() + dx, (*ite)->y() + dy, (*ite)->z() + dz);
	}
	/** 修改几何顶点坐标 **/
	for (auto v : f->vertices()) {
		for (std::size_t j = 0; j < v->getNumMeshVertices(); ++j) {
			GPoint p(v->x() + dx, v->y() + dy, v->z() + dz);
			v->setPosition(p);
		}
	}
	/** 更新distance field **/
	updateField();
}


/*!
 * @brief      移动某个区域内的网格。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  GRegion * r  
 * @param[in]  double dx  
 * @param[in]  double dy  
 * @param[in]  double dz  
 * @return     void  
**/
void Relay::moveRegion(GRegion* r, double dx, double dy, double dz) {
	//    printf("r->getNumMeshVertices():%llu\n",r->getNumMeshVertices());
	//    printf("r->embeddedVertices():%llu\n",r->embeddedVertices().size());
	//    printf("r->vertices():%llu\n",r->vertices().size());

    /** 移动内部的节点 **/
	std::size_t num = r->getNumMeshVertices();
	for (std::size_t j = 0; j < num; ++j) {
		r->getMeshVertex(j)->setXYZ(r->getMeshVertex(j)->x() + dx,
			r->getMeshVertex(j)->y() + dy,
			r->getMeshVertex(j)->z() + dz);
	}

	/** 修改边界面上分网节点坐标 **/
	std::set<MVertex *, MVertexLessThanNum> all_vertices;
	for (auto f : r->faces()) {
		/** 修改内部分网节点坐标 **/
		num = f->getNumMeshVertices();
		for (std::size_t j = 0; j < num; ++j) {
			f->getMeshVertex(j)->setXYZ(f->getMeshVertex(j)->x() + dx,
				f->getMeshVertex(j)->y() + dy,
				f->getMeshVertex(j)->z() + dz);
		}
		/** 由于两个面相邻的话，分别移动边界上的点会造成重复移动 **/
		for (auto e : f->edges()) {
			for (auto line : e->lines) {
				MVertex *v1 = line->getVertex(0);
				MVertex *v2 = line->getVertex(1);

				all_vertices.insert(v1);
				all_vertices.insert(v2);
			}
		}
	}
	/** 移动面的边上的分网点**/
	for (std::set<MVertex *, MVertexLessThanNum>::iterator ite = all_vertices.begin(); ite != all_vertices.end(); ++ite) {
		(*ite)->setXYZ((*ite)->x() + dx, (*ite)->y() + dy, (*ite)->z() + dz);
	}
	/** 修改几何顶点坐标 **/
	for (auto v : r->vertices()) {
		for (std::size_t j = 0; j < v->getNumMeshVertices(); ++j) {
			GPoint p(v->x() + dx, v->y() + dy, v->z() + dz);
			v->setPosition(p);
		}
	}
	/** 更新distance field **/
	updateField();
}

/**
 * @brief 旋转某个区域内的网格。正方向为右手握住轴的方向。
 *
 * @param r
 * @param dangle -180~180
 * @param x1 向量的位置
 * @param y1
 * @param z1
 * @param v1 旋转单位向量，过原点。
 * @param v2
 * @param v3
 */
void Relay::rotateRegion(GRegion* r, double dangle, double x1, double y1, double z1, double v1, double v2, double v3) {
	Point3f p;
	/** 移动内部的节点 **/
	std::size_t num = r->getNumMeshVertices();
	for (std::size_t j = 0; j < num; ++j) {
		p = RotateByVector(dangle, r->getMeshVertex(j)->x() - x1, r->getMeshVertex(j)->y() - y1, r->getMeshVertex(j)->z() - z1, v1, v2, v3);
		r->getMeshVertex(j)->setXYZ(p.x + x1, p.y + y1, p.z + z1);
	}

	/** 修改边界面上分网节点坐标 **/
	std::set<MVertex *, MVertexLessThanNum> all_vertices;
	for (auto f : r->faces()) {
		/** 修改内部分网节点坐标 **/
		num = f->getNumMeshVertices();
		for (std::size_t j = 0; j < num; ++j) {
			p = RotateByVector(dangle, f->getMeshVertex(j)->x() - x1, f->getMeshVertex(j)->y() - y1, f->getMeshVertex(j)->z() - z1, v1, v2, v3);
			f->getMeshVertex(j)->setXYZ(p.x + x1, p.y + y1, p.z + z1);
		}
		/** 由于两个面相邻的话，分别移动边界上的点会造成重复移动 **/
		for (auto e : f->edges()) {
			for (auto line : e->lines) {
				MVertex *v1 = line->getVertex(0);
				MVertex *v2 = line->getVertex(1);

				all_vertices.insert(v1);
				all_vertices.insert(v2);
			}
		}
	}
    /** 移动面的边上的分网点 **/
	for (std::set<MVertex *, MVertexLessThanNum>::iterator ite = all_vertices.begin(); ite != all_vertices.end(); ite++) {
		p = RotateByVector(dangle, (*ite)->x() - x1, (*ite)->y() - y1, (*ite)->z() - z1, v1, v2, v3);
		(*ite)->setXYZ(p.x + x1, p.y + y1, p.z + z1);
	}
	/** 修改几何顶点坐标 **/
	for (auto v : r->vertices()) {
		num = v->getNumMeshVertices();
		for (std::size_t j = 0; j < num; j++) {
			p = RotateByVector(dangle, v->x() - x1, v->y() - y1, v->z() - z1, v1, v2, v3);
			GPoint p1(p.x + x1, p.y + y1, p.z + z1);
			v->setPosition(p1);
		}
	}
	/** 更新distance field **/
	updateField();
}

/**
 * @brief 空间点绕任意轴的旋转公式。
 *
 * @param theta 角度
 * @param old_x
 * @param old_y
 * @param old_z
 * @param vx 单位向量。
 * @param vy
 * @param vz
 * @return Point3f
 */
Point3f Relay::RotateByVector(double theta, double old_x, double old_y, double old_z, double vx, double vy, double vz) {
	double r = theta * M_PI / 180;
	double c = cos(r);
	double s = sin(r);

	double new_x = (vx*vx*(1 - c) + c) * old_x + (vx*vy*(1 - c) - vz * s) * old_y + (vx*vz*(1 - c) + vy * s) * old_z;
	double new_y = (vy*vx*(1 - c) + vz * s) * old_x + (vy*vy*(1 - c) + c) * old_y + (vy*vz*(1 - c) - vx * s) * old_z;
	double new_z = (vx*vz*(1 - c) - vy * s) * old_x + (vy*vz*(1 - c) + vx * s) * old_y + (vz*vz*(1 - c) + c) * old_z;
	return Point3f(new_x, new_y, new_z);
}


/*!
 * @brief      对模型进行重分网，对衔铁进行移动，对可压缩区域进行重分网。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  double dx  
 * @param[in]  double dy  
 * @return     void  
**/
void Relay::remesh(double dx, double dy) {
    printf("mesh step %d\n", current_step);
    char nameMsh[256];
	sprintf(nameMsh, "%s_%02d.msh", fileName, current_step);
	/** 未发生位移就不要分了 **/
	if (fabs(dx) < 1e-10 && fabs(dy) < 1e-10) {
		printf("Skipping remeshing.\n");
		gmsh::write(nameMsh);
		return;
	}

	/** 删掉空气的分网 **/
	GFace* f_xiantie = nullptr;
	GFace* f_air = nullptr;
	for (GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it) {
		if ((*it)->tag() == tag_remesh) {
			f_air = (*it);
		}
		if ((*it)->tag() == tag_xiantie) {
			f_xiantie = (*it);
		}
	}    

	/** 移动衔铁的分网 **/
	if (f_xiantie) {
        printf("moving armature region %d mesh dx=%lf,dy=%lf...\n", tag_xiantie, dx, dy);
		moveFace(f_xiantie, dx, dy, 0);
	}
	
	/** 判断形变网格的质量 **/
	RemeshOrMorph = true;

	if (f_air && RemeshOrMorph) {
		/** 删除空气的分网 **/
        printf("deleting air surface %d\n", tag_remesh);
		deleteFaceMesh(f_air);
		/** 对空气进行重分网 **/
        printf("remesh air domain...\n");
		f_air->mesh(true);
	}

	gmsh::write(nameMsh);

    printf("Finish remesh 2d\n");
}


/*!
 * @brief      对3D分网进行平行移动。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  double dx  移动的距离，默认大小与坐标轴一致。
 * @param[in]  double dy  
 * @param[in]  double dz  
 * @return     void  
**/
void Relay::remesh3DParallel(double dx, double dy, double dz) {
    printf("mesh step %d\n", current_step);
    char nameMsh[256];
	sprintf(nameMsh, "%s_%02d.msh", fileName, current_step);
	/** 未发生位移就不要分了 **/
	if (dx == 0 && dy == 0 && dz == 0) {
        printf("Skipping remeshing.\n");
        gmsh::write(nameMsh);
		return;
	}

	GRegion* r_xiantie = nullptr;
	GRegion* r_air = nullptr;
    for (GModel::riter it = model->firstRegion (); it != model->lastRegion(); ++it) {
		//printf("rigion tag:%d\n", (*it)->tag());
		if ((*it)->tag() == tag_remesh) {
			r_air = (*it);
		}
		if ((*it)->tag() == tag_xiantie) {
			r_xiantie = (*it);
		}
	}    
	
	/** 移动衔铁的分网 **/
	if (r_xiantie) {
        printf("moving armature region mesh %d dx=%10.8e, dy=%10.8e, dz=%10.8e...\n", tag_xiantie, dx, dy, dz);
		moveRegion(r_xiantie, dx, dy, dz);
	}

	/** 判断形变网格的质量 **/
    RemeshOrMorph = true;

	if (r_air && RemeshOrMorph) {
		/** 删除空气的分网 **/
        printf("deleting air region %d\n", tag_remesh);
        deleteRegionMesh(r_air);
		/** 对空气进行重分网 **/
        printf("remesh air domain %d\n", tag_remesh);
        /** gmsh并没有三维的分网接口 **/
        meshGRegionNetgen(r_air);
        model->setAllVolumesPositive();
	}

	gmsh::write(nameMsh);

    printf("Finish remesh3DParallel\n");
}


/*!
 * @brief      对分网进行旋转移动。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  double dangle  移动的角度，方向是右手握住轴的方向，即逆时针为正。
 * @param[in]  double x1  轴上的点x
 * @param[in]  double y1  轴上的点y
 * @param[in]  double z1  轴上的点z
 * @param[in]  double v1  轴向量x
 * @param[in]  double v2  轴向量y
 * @param[in]  double v3  轴向量z
 * @return     void  
**/
void Relay::remesh3DRotate(double dangle, double x1, double y1, double z1, double v1, double v2, double v3) {
    printf("mesh step %d\n", current_step);
    char nameMsh[256];
	sprintf(nameMsh, "%s_%02d.msh", fileName, current_step);
	/** 未发生位移就不要分了 **/
	if (dangle == 0) {
        printf("Skipping remeshing.\n");
		gmsh::write(nameMsh);
		return;
	}

	GRegion* r_xiantie = nullptr;
	GRegion* r_air = nullptr;
	for (GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it) {
		//printf("region tag:%d\n", (*it)->tag());
		if ((*it)->tag() == tag_remesh) {
			r_air = (*it);
		}
		if ((*it)->tag() == tag_xiantie) {
			r_xiantie = (*it);
		}
	}    
	
	/** 移动衔铁的分网 **/
	if (r_xiantie) {
        printf("rotate armature region mesh %d dangle=%lf...\n", tag_xiantie, dangle);
		rotateRegion(r_xiantie, dangle, x1, y1, z1, v1, v2, v3);
	}
	
	/** 判断形变网格的质量 **/
    RemeshOrMorph = false;

	if (r_air && RemeshOrMorph) {
		/** 删除空气的分网 **/
        printf("deleting air region %d\n", tag_remesh);
        deleteRegionMesh(r_air);
		/** 对空气进行重分网 **/
        printf("remesh air domain %d\n", tag_remesh);
        /** gmsh并没有三维的分网接口 **/
        meshGRegionNetgen(r_air);
        model->setAllVolumesPositive();
	}	

	gmsh::write(nameMsh);

    printf("Finish remesh3DRotate\n");
}


/*!
 * @brief      更新gmsh的分网大小控制。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     void  
**/
void Relay::updateField() {
	/** 更新distance field **/
	FieldManager* fM = GModel::current()->getFields();
	std::map<int, Field *>::iterator iter;
	for (iter = fM->begin(); iter != fM->end(); ++iter) {
		iter->second->update_needed = true;
		iter->second->update();
	}
}


/*!
 * @brief      检查生成的连通矩阵是否正确
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     bool  
**/
bool Relay::checkConnectivityMatrix() {
	int edge_size = allEdges.size();
	std::vector<double> allEdgesLength(edge_size);

	double x1, y1, z1, x2, y2, z2;
	int start, end;
	/** 计算棱长 **/
	for (int i = 0; i < edge_size; ++i) {
		start = allEdges[i].start;
		end = allEdges[i].end;
		x1 = pmeshnode[start].x;
		y1 = pmeshnode[start].y;
		z1 = pmeshnode[start].z;
		x2 = pmeshnode[end].x;
		y2 = pmeshnode[end].y;
		z2 = pmeshnode[end].z;
		allEdgesLength[i] = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
	}

	/** 根据原始信息计算棱长 **/
	int n[10];
    for (int i = 0; i < num_tet; ++i) {
		int i1 = num_ele - num_tet + i;
		CElement* ele = pmeshele + i1;

		for (int in = 0; in < 10; ++in) {
			n[in] = ele->n[in];
		}
		double l = 0;
		for (int in = 4; in < 10; ++in) {
			start = n[tet_n[in - 4][0]];
			end = n[tet_n[in - 4][1]];
			x1 = pmeshnode[start].x;
			y1 = pmeshnode[start].y;
			z1 = pmeshnode[start].z;
			x2 = pmeshnode[end].x;
			y2 = pmeshnode[end].y;
			z2 = pmeshnode[end].z;
			l = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
			if (std::abs(l - allEdgesLength[n[in]]) > 1e-10) {
				printf("%lf\t", std::abs(l - allEdgesLength[n[in]]));
				return false;
			}
		}
	}
	return true;
}


/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     void  
**/
void Relay::run() {

}


/*!
 * @brief      初始化求解中的变量。
 * @author     Poofee
 * @date       2020年08月30日
 * @param[out] 
 * @param[in]  int size  
 * @return     void  
**/
void Relay::initPhysicalVariables(int size) {
	totaltimes.resize(size + 1);
	totaltimes[0] = 0;
	Ddisplacements.resize(size + 1);
	Ddisplacements[0] = 0;
	displacements.resize(size + 1);
	displacements[0] = 0;
	velocities.resize(size + 1);
	velocities[0] = 0;
	accelerations.resize(size + 1);
	accelerations[0] = 0;
	PhiCoil.resize(size + 1);
	PhiCoil[0] = 0;
	ICoil.resize(size + 1);
	ICoil[0] = 0;
	magforces.resize(size + 1);
	magforces[0] = 0;
	UCoil.resize(size + 1);
    UCoil[0] = 0;
    IterationStep.resize(size + 1);
    SimulationTime.resize(size + 1);
}

/*!
 * @brief      更新稀疏矩阵位置向量。
 * @author     Poofee
 * @date       2020年11月14日
 * @param[out]
 * @param[in]  int pos
 * @param[in]  int row
 * @param[in]  int col
 * @param[in]  double vals
 * @return     void
**/
void Relay::insertSparseMatrixData(int pos, int row, int col, double values) {
    locs(0, pos) = row;
    locs(1, pos) = col;
    vals(0, pos) = values;
}


/*!
 * @brief      给区域index设置材料。
 * @author     Poofee
 * @date       2020年09月18日
 * @param[out] 
 * @param[in]  int index  
 * @param[in]  CMaterial * mat  
 * @return     void  
**/
void Relay::setMaterial(int index, CMaterial* mat) {
	materialMap[index] = mat;
}


/*!
 * @brief      获取index区域的材料。
 * @author     Poofee
 * @date       2020年09月18日
 * @param[out] 
 * @param[in]  int domainIndex  
 * @return     CMaterial*  
**/
CMaterial* Relay::getMaterial(int domainIndex) {
	return materialMap[domainIndex];
}

/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @return     void  
**/
void Relay::stepIncrement() {
    current_step++;
}
/*!
 * @brief      设置prefix名。
 * @author     Poofee
 * @date       2020年11月5日
 * @param[out]
 * @param[in]  char * fn
 * @return     void
**/
void Relay::setPrefixName(char* fn) {
    sprintf(prefixName, "%s", fn);
}

/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int index  
 * @return     double  
**/
double Relay::calcMagForce(int index, FORCEMETHOD method) {
    printf("Must be implemented in subclass.\n");
	return 0;
}

/*!
 * @brief      
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  int index  
 * @return     double  
**/
double Relay::calcMagTorque(int index, FORCEMETHOD method) {
    printf("Must be implemented in subclass.\n");
	return 0;
}


/*!
 * @brief      输出结果数据。MATLAB格式。
 * @author     Poofee
 * @date       2020年10月3日
 * @param[out] 
 * @return     void  
**/
void Relay::outputResults() {
	/** 保存到txt **/
	char fn[256];
	sprintf(fn, "%sresults.m", fileName);
	FILE* fp;
	fp = fopen(fn, "w");
	fprintf(fp, "%%output by FEEM\n");
	fprintf(fp, "%%timesteps, displacements, Ddisplacements, velocities, accelerations,PhiCoil, ICoil, magforces, UCoil\n");

	fprintf(fp, "results = [\n"); 
	for (int i = 0; i < timesteps.size(); ++i) {
		fprintf(fp, "%10.8e,", totaltimes[i]);
		fprintf(fp, "%10.8e,", displacements[i]);
		fprintf(fp, "%10.8e,", velocities[i]);
		fprintf(fp, "%10.8e,", accelerations[i]);
		fprintf(fp, "%10.8e,", magforces[i]);
		fprintf(fp, "%10.8e,", ICoil[i]);
		fprintf(fp, "%10.8e,", PhiCoil[i]);
		fprintf(fp, "%10.8e,", UCoil[i]);
		fprintf(fp, "%10.8e,", Ddisplacements[i]);
        fprintf(fp, "%10.8e,", IterationStep[i]);
        fprintf(fp, "%10.8e,", SimulationTime[i]);
		fprintf(fp, "; \n");
	}
	fprintf(fp, "];\n");

	fprintf(fp, "subplot(2,3,1);hold on;\n");
	fprintf(fp, "plot(results(:,1),results(:,2),'*-');\n");
	//fprintf(fp, "plot(results(:,1),results(:,3),'*-');\n");
	fprintf(fp, "title(\"%s\");\n\n","displacement");

	fprintf(fp, "subplot(2,3,2);hold on;\n");
	fprintf(fp, "plot(results(:,1),results(:,3),'*-');\n");
	fprintf(fp, "title(\"%s\");\n\n", "velocities");

	fprintf(fp, "subplot(2,3,3);hold on;\n");
	fprintf(fp, "plot(results(:,1),results(:,4),'*-');\n");
	fprintf(fp, "title(\"%s\");\n\n", "accelerations");

	fprintf(fp, "subplot(2,3,4);hold on;\n");
	fprintf(fp, "plot(results(:,1),results(:,5),'*-');\n");
	fprintf(fp, "title(\"%s\");\n\n", "magforces");

	fprintf(fp, "subplot(2,3,5);hold on;\n");
	fprintf(fp, "plot(results(:,1),results(:,6),'*-');\n");
	//fprintf(fp, "plot(results(:,1),results(:,9),'*-');\n");
	fprintf(fp, "title(\"%s\");\n\n", "ICoil");

	fprintf(fp, "subplot(2,3,6);hold on;\n");
	fprintf(fp, "plot(results(:,1),results(:,7),'*-');\n");
	fprintf(fp, "title(\"%s\");\n", "PhiCoil");

	fclose(fp);
}

/*!
 * @brief      写vtk文件
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  char * fname  
 * @param[in]  int num_cell  
 * @param[in]  int ele_type  
 * @param[in]  int ele_size  
 * @return     void  
**/
void Relay::writeVtkFile(char* fname, int num_cell, int ele_type, int ele_size) {
	FILE*  fp;
	fp = fopen(fname, "w");

	/*
	 1: points
	 3: line
	 5: Triangular element
	 9: Quadrilateral element
	10: Tetrahedral element
	12: Hexahedral element
	13: Triangular prism element
	14: Pyramid element
	*/
	/** 数据版本声明 **/
	fprintf(fp, "# vtk DataFile Version 2.0\n");
	/** 标题 **/
	fprintf(fp, "vtk title\n");
	/** 文件格式声明 **/
	fprintf(fp, "ASCII\n");
	/** 几何拓扑结构 **/
	fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
	/** 节点 **/
	fprintf(fp, "\nPOINTS %d double\n", num_pts);
	for (int i = 0; i < num_pts; ++i) {
		fprintf(fp, "%lf %lf %lf\n", pmeshnode[i].x, pmeshnode[i].y, pmeshnode[i].z);
	}
	/** 单元 **/
	fprintf(fp, "\nCELLS %d %d\n", num_cell, (ele_size + 1)*num_cell);
	for (int i_tri = 0; i_tri < num_cell; ++i_tri) {
		int i1 = num_ele - num_cell + i_tri;
		fprintf(fp, "%d", ele_size);
		for (int i = 0; i < ele_size; ++i) {
			fprintf(fp, " %d", pmeshele[i1].n[i]);
		}
		fprintf(fp, "\n");
	}
	/** 单元类型 **/
	fprintf(fp, "\nCELL_TYPES %d\n", num_cell);
	for (int i = 0; i < num_cell; ++i) {
		fprintf(fp, "%d\n", ele_type);
	}

	fclose(fp);
}


/*!
 * @brief      写出边界。
 * @author     Poofee
 * @date       2020年10月5日
 * @param[out] 
 * @param[in]  char fn[]  
 * @return     void  
**/
void Relay::testBoundaryResults(char fn[]) {
	/** 调试用，验证结果是否正确 **/
	FILE*  fp;
	fp = fopen(fn, "w");

	/*
		 1: points
		 3: line
		 5: Triangular element
		 9: Quadrilateral element
		10: Tetrahedral element
		12: Hexahedral element
		13: Triangular prism element
		14: Pyramid element
	*/
	/** 数据版本声明 **/
	fprintf(fp, "# vtk DataFile Version 2.0\n");
	/** 标题 **/
	fprintf(fp, "vtk title\n");
	/** 文件格式声明 **/
	fprintf(fp, "ASCII\n");
	/** 几何拓扑结构 **/
	fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
	/** 节点 **/
	//fprintf(fp, "\nPOINTS %d float\n", boundaryPoints.size());
	//for (int i = 0; i < boundaryPoints.size(); ++i) {
	//	fprintf(fp, "%lf %lf %lf\n", pmeshnode[boundaryPoints[i]].x, pmeshnode[boundaryPoints[i]].y, pmeshnode[boundaryPoints[i]].z);
	//}
	fprintf(fp, "\nPOINTS %d float\n", num_pts);
	for (int i = 0; i < num_pts; ++i) {
		fprintf(fp, "%lf %lf %lf\n", pmeshnode[i].x, pmeshnode[i].y, pmeshnode[i].z);
	}
	/** 单元 **/
	int num_cell = boundaryEdges.size();
	int ele_size = 1;
	int ele_type = 1;
	fprintf(fp, "\nCELLS %d %d\n", num_cell + boundaryPoints.size(), (2 + 1)*num_cell + boundaryPoints.size() * 2);
	for (int i_tri = 0; i_tri < boundaryPoints.size(); ++i_tri) {
		fprintf(fp, "%d", ele_size);
		fprintf(fp, " %d", boundaryPoints[i_tri]);
		fprintf(fp, "\n");
	}
	ele_size = 2;
	ele_type = 3;
	for (int i_tri = 0; i_tri < num_cell; ++i_tri) {
		fprintf(fp, "%d", ele_size);
		fprintf(fp, " %d", boundaryEdges[i_tri].start);
		fprintf(fp, " %d", boundaryEdges[i_tri].end);
		fprintf(fp, "\n");
	}

	/** 单元类型 **/
	fprintf(fp, "\nCELL_TYPES %d\n", num_cell + boundaryPoints.size());
	ele_type = 1;
	for (int i_tri = 0; i_tri < boundaryPoints.size(); ++i_tri) {
		fprintf(fp, "%d\n", ele_type);
	}
	ele_type = 3;
	for (int i = 0; i < num_cell; ++i) {
		fprintf(fp, "%d\n", ele_type);
	}
    fclose(fp);
}

/*!
 * @brief      打印当前时间。
 * @author     Poofee
 * @date       2020年11月09日
 * @param[out]
 * @return     void
**/
void Relay::printCurrentTime() {
    time_t timep;

    time(&timep);

    struct tm *t = NULL;

    t = localtime(&timep);

    printf("\n%04d-%02d-%02d %02d:%02d:%02d\n", 1900 + t->tm_year, 1 + t->tm_mon, t->tm_mday, t->tm_hour, t->tm_min, t->tm_sec);
}

/*!
 * @brief      静态参数初始化
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay::StaticInitial() {

}

/*!
 * @brief      瞬态参数初始化
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay::TransientInitial() {

}

/*!
 * @brief      静磁场分析
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay::StaticAnalysis() {
	printf("Please finish it in subclass\n");
}

/*!
 * @brief      瞬态磁场分析
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay::TransientAnalysis() {
	printf("Please finish it in subclass\n");
}

